PSSFSS Examples

Symmetric Strip Grating

This example consists of a symmetric strip grating, i.e. a grating where the strip width is half the unit cell period $P$:

Only three of the infinite number of strips in the grating are shown, and they extend infinitely to the left and right. The grating lies in the $z=0$ plane with free space on both sides. The shaded areas represent metalization. The dashed lines show two possible choices for the unit cell location: "J" for a formulation in terms of electric surface currents, and "M" for magnetic surface currents.

For normal incidence there is a closed-form solution due to Weinstein, but for a more recent reference one can find the solution in Problem 10.6 of R. E. Collin, Field Theory of Guided Waves, Second Ed., IEEE Press, 1991. Here is the code for computing the exact solution:

"""
    grating(kP, nterms=30) -> (Γ, T)

Compute the normal incidence reflecton and transmission coefficients of a symmetric grid of
zero-thickness conducting strips.  The product of the period of the strips and the incident
electric field wavenumber is `kP` (dimensionless).  The incident electric field is
perpendicular to the direction along the axis of the strips.  The series have been
accelerated by applying a Kummer's transformation, using the first two terms in the Maclaurin
series for the inverse sin function.  `kP` must be in the half-open interval [0,1). The
default number of summed terms `nterms` yields better than 10 digits of accuracy over the
interval [0.01,0.99].
"""
function grating(kP; nterms=30)
    sum1 = 1.3862943611198906 # \sum_{n=1}^{\infty} 1/(n-1/2) - 1/n = log(4)
    sum3 = 7.2123414189575710 # \sum_{n=1}^{\infty} (n-1/2)^{-3} - n^{-3} = 6 * \zeta(3)
    x = kP/(4π)
    θ = x*sum1 + x^3/6 * sum3
    for n = 1:nterms
        xonmhalf = x/(n - 0.5)
        xon = x/n
        term = asin(xonmhalf) - (xonmhalf + (xonmhalf)^3/6) -
              (asin(xon) - (xon + xon^3/6))
        θ += term
    end
    Γ = sin(θ) * cis(-π/2 - θ)
    T = 1 + Γ
    return (Γ, T)
end

Main.ex-symmetric_strip.grating

Note that using the extension of Babinet's Principle for electromagnetic fields this also provides the solution (upon appropriate interchange and sign change of the coefficients) for the case where the incident wave polarization is parallel to the direction of the strips.

Here is the PSSFSS code to analyze this structure using electric currents as the unknowns. We will scale the geometry so that the frequency in GHz is numerically equal to the period of the strips measured in wavelengths.

using Plots, PSSFSS
c = 11.802852677165355 # light speed [inch*GHz]
period = c  # so the period/wavelength = freq in GHz
Py = period
Ly = period/2
Px = Lx = Ly/10 # Infinite in x direction so this can be anything
Ny = 60
Nx = round(Int, Ny*Lx/Ly)
sheet = rectstrip(;Px, Py, Lx, Ly, Nx, Ny, units=inch)
flist = 0.02:0.02:0.98
steering = (θ=0, ϕ=0)
strata = [Layer()
          sheet
          Layer()]
results_j = analyze(strata, flist, steering, showprogress=false,
                    resultfile=devnull, logfile=devnull);
p1 = plot(sheet)
p2 = plot(sheet, unitcell=true)
ptitle = plot(title = "Symmetric Strip Triangulation",
             grid = false, showaxis = false, xtick=[], ytick=[],
             bottom_margin = -50Plots.px)
plot(ptitle, p1, p2, layout = @layout([A{0.09h}; [B C]]))

Note that setting Lx = Px causes the strip to fully occupy the x-extent of the unit cell. PSSFSS automatically ensures that the triangle edges at these unit cell boundaries define basis functions that satisfy the Floquet (phase shift) boundary conditions, so that currents are free to flow across these unit cell boundaries.

We can also analyze the same structure using magnetic currents in the areas free of metalization as the unknowns:

sheet = rectstrip(;class='M', Px, Py, Lx, Ly, Nx, Ny, units=inch)
strata = [Layer()
          sheet
          Layer()]
results_m = analyze(strata, flist, steering, showprogress=false,
                    resultfile=devnull, logfile=devnull);

Each 50-frequency run of analyze takes about 14 seconds for this geometry of 720 triangles on my machine. More detailed timing information is available in the log file (which is omitted for generating this documentation).

We will compare the PSSFSS results to the analytic solution:

# Generate exact results:
rt = grating.(2π*flist)
rperp_exact = first.(rt)
tperp_exact = last.(rt)
rpar_exact = -tperp_exact
tpar_exact = -rperp_exact;

Obtain PSSFSS results for electric and magnetic currents:

outrequest = @outputs s11(v,v) s21(v,v) s11(h,h) s21(h,h)
rperp_j, tperp_j, rpar_j, tpar_j =
      collect.(eachcol(extract_result(results_j, outrequest)))
rperp_m, tperp_m, rpar_m, tpar_m =
      collect.(eachcol(extract_result(results_m, outrequest)));

Generate the comparison plots:

angdeg(z) = rad2deg(angle(z)) # Convenience function

p1 = plot(title = "Perpendicular Reflection Magnitude",
          xlabel = "Period (wavelengths)",
          ylabel = "Coefficient Magnitude",
          legend=:topleft)
plot!(p1, flist, abs.(rperp_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(rperp_j), label="PSSFSS J")
plot!(p1, flist, abs.(rperp_m), label="PSSFSS M")

p2 = plot(title = "Perpendicular Reflection Phase",
          xlabel = "Period (wavelengths)",
          ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(rperp_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(rperp_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(rperp_m), label="PSSFSS M")

p1 = plot(title = "Parallel Reflection Magnitude",
          xlabel = "Period (wavelengths)",
          ylabel = "Coefficient Magnitude")
plot!(p1, flist, abs.(rpar_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(rpar_j), label="PSSFSS J")
plot!(p1, flist, abs.(rpar_m), label="PSSFSS M")

p2 = plot(title = "Parallel Reflection Phase",
          xlabel = "Period (wavelengths)",
          ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(rpar_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(rpar_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(rpar_m), label="PSSFSS M")

Now look at the transmission coefficients:

p1 = plot(title = "Perpendicular Transmission Magnitude",
          xlabel = "Period (wavelengths)",
          ylabel = "Coefficient Magnitude")
plot!(p1, flist, abs.(tperp_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(tperp_j), label="PSSFSS J")
plot!(p1, flist, abs.(tperp_m), label="PSSFSS M")

p2 = plot(title = "Perpendicular Transmission Phase",
          xlabel = "Period (wavelengths)",
          ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(tperp_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(tperp_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(tperp_m), label="PSSFSS M")

p1 = plot(title = "Parallel Transmission Magnitude",
          xlabel = "Period (wavelengths)",
          ylabel = "Coefficient Magnitude", legend=:topleft)
plot!(p1, flist, abs.(tpar_exact), ls=:dash, label="Exact")
plot!(p1, flist, abs.(tpar_j), label="PSSFSS J")
plot!(p1, flist, abs.(tpar_m), label="PSSFSS M")

p2 = plot(title = "Parallel Transmission Phase",
          xlabel = "Period (wavelengths)",
          ylabel = "Phase (deg)")
plot!(p2, flist, angdeg.(tpar_exact), ls=:dash, label="Exact")
plot!(p2, flist, angdeg.(tpar_j), label="PSSFSS J")
plot!(p2, flist, angdeg.(tpar_m), label="PSSFSS M")

Conclusion

Although good agreement is obtained, as expected the best agreement between all three results occurs for the lowest frequencies, where the triangles are smallest in terms of wavelength. This emphasizes the fact that it is necessary for the user to check that enough triangles have been requested for good convergence over the frequency band of interest. This example is an extremely demanding case in terms of bandwidth, as the ratio of maximum to minimum frequency here is $0.98/0.02 = 49:1$

Resistive Square Patch

This example will demonstrate the ability of PSSFSS to accurately model finite conductivity of FSS metalization. It consists of a square finitely conducting patch in a square lattice. It is taken from a paper by Alon S. Barlevy and Yahya Rahmat-Samii, "Fundamental Constraints on the Electrical Characteristics of Frequency Selective Surfaces", Electromagnetics, vol. 17, 1997, pp. 41-68. This particular example is from Section 3.2, Figures 7 and 8. We will compare PSSFSS results to those digitized from the cited figures.

We start by defining a function that creates a patch of the desired sheet resistance:

using Plots, PSSFSS
patch(R) = rectstrip(Nx=10, Ny=10, Px=1, Py=1, Lx=0.5, Ly=0.5, units=cm, Rsheet=R)
plot(patch(0), unitcell=true)

The patches measure 0.5 cm on a side and lie in a square lattice of period 1 cm. Now we perform the analysis, looping over the desired values of sheet resistance.

steering = (ϕ=0, θ=0)
flist = 1:0.5:60
Rs = [0, 10, 30, 100]
calculated = zeros(length(flist), length(Rs)) # preallocate storage
outputs = @outputs s11mag(v,v)
for (i,R) in pairs(Rs)
    strata = [Layer(), patch(R), Layer()]
    results = analyze(strata, flist, steering, showprogress=false,
                      logfile=devnull, resultfile=devnull)
    calculated[:,i] = extract_result(results, outputs)
end

Looping over the four sheet resistance values, each evaluated at 119 frequencies required approximately 20 seconds on my machine.

We plot the results, including those digitized from the paper for comparison:

using DelimitedFiles
markers = (:diamond, :utriangle, :square, :xcross)
colors = (:blue, :red, :green, :black)
p = plot(xlim=(-0.01,60.01), xtick = 0:10:60, ylim=(-0.01,1.01), ytick=0:0.1:1,
         xlabel="Frequency (GHz)", ylabel="Reflection Coefficient Magnitude",
         title = "Resistive Square Patch",
         legend=:topright)
for (i,R) in pairs(Rs)
    scatter!(p, flist, calculated[:,i], label="PSSFSS $R Ω", ms=2, shape=markers[i], color=colors[i])
    data = readdlm("../src/assets/barlevy_patch_digitized_$(R)ohm.csv", ',')
    plot!(p, data[:,1], data[:,2], label="Barlevy $R Ω", ls=:dash, color=colors[i])
end
p

Conclusion

PSSFSS results are indistinguishable from those reported in the cited paper.

Cross on Dielectric Substrate

This example is also taken from the paper by Alon S. Barlevy and Yahya Rahmat-Samii, "Fundamental Constraints on the Electrical Characteristics of Frequency Selective Surfaces", Electromagnetics, vol. 17, 1997, pp. 41-68. This particular example is from Section 3.2, Figures 7 and 8. It also appeared at higher resolution in Barlevy's PhD dissertation from which the comparison curves were digitized.

We use the loadedcross element where we choose w > L2/2, so that the Cross is "unloaded", i.e. the center section is filled in with metalization:

using Plots, PSSFSS, DelimitedFiles
sheet = loadedcross(w=1.0, L1=0.6875, L2=0.0625, s1=[1.0,0.0],
                    s2=[0.0,1.0], ntri=600, units=cm)
plot(sheet, unitcell=true)

A few things to note. First, the mesh is unstructured. So there are no redundant triangle face-pairs that PSSFSS can exploit to reduce execution time. Second, the number of triangle faces generated is only approximately equal to the requested value of 600. This can be verified by entering the Julia variable sheet at the REPL (i.e. the Julia prompt):

sheet

RWGSheet: style=loadedcross, class=J, 404 nodes, 1054 edges, 651 faces

The cross FSS is etched on a dielectric sheet of thickness 3 mm. The dielectric constant is varied over the values 1, 2, and 4 to observe the effect on the resonant frequency. Following the reference, the list of analysis frequencies is varied slightly depending on the value of dielectric constant:

resultsstack = Any[]
steering = (ϕ=0, θ=0)
for eps in [1, 2, 4]
    strata = [  Layer()
                sheet
                Layer(ϵᵣ=eps, width=3mm)
                Layer()
             ]
    if eps == 1
        flist = 1:0.2:30
    elseif eps == 2
        flist = 1:0.2:26
    else
        flist = 1:0.2:20
    end
    results = analyze(strata, flist, steering, showprogress=false,
                      resultfile=devnull, logfile=devnull)
    push!(resultsstack, results)
end

The above loop requires about 80 seconds of execution time on my machine. Compare PSSFSS results to those digitized from the dissertation figure:

col=[:red,:blue,:green]
p = plot(xlim=(0.,30), xtick = 0:5:30, ylim=(0,1), ytick=0:0.1:1,
         xlabel="Frequency (GHz)", ylabel="Reflection Coefficient Magnitude",
         legend=:topleft, lw=2)
for (i,eps) in enumerate([1,2,4])
    data = extract_result(resultsstack[i],  @outputs FGHz s11mag(v,v))
    plot!(p, data[:,1], data[:,2], label="PSSFSS ϵᵣ = $eps", lc=col[i])
    data = readdlm("../src/assets/barlevy_diss_eps$(eps).csv", ',')
    plot!(p, data[:,1], data[:,2], label="Barlevy ϵᵣ = $eps", lc=col[i], ls=:dot)
end
p

Conclusion

PSSFSS results agree very well with those of the cited reference, especially when accounting for the fact that the reference results are obtained by digitizing a scanned figure.

Square Loop Absorber

This example is from Figure 7 of Costa and Monorchio: "A frequency selective radome with wideband absorbing properties", IEEE Trans. AP-S, Vol. 60, no. 6, June 2012, pp. 2740–2747. It shows how one can use the polyring function to model square loop elements. Three different designs are examined that employ different loop thicknesses and different values of sheet resistance. We compare the reflection coefficient magnitudes computed by PSSFSS with those digitized from the cited figure when the sheet is suspended 5 mm above a ground plane, hence we will also make use of the pecsheet function.

using Plots, PSSFSS, DelimitedFiles
D = 11 # Period of square lattice (mm)
r_outer = √2/2 * D/8 * [5,6,7] # radii of square outer vertices
thickness = D/16 * [1,2,3]
r_inner = r_outer - √2 * thickness  # radii of square inner vertices
Rs = [15,40,70] # Sheet resistance (Ω/□)
labels = ["Thin", "Medium", "Thick"]
colors = [:green, :blue, :red]
p = plot(title="Costa Absorber", xlim=(0,25),ylim=(-35,0),xtick=0:5:25,ytick=-35:5:0,
         xlabel="Frequency (GHz)", ylabel="Reflection Magnitude (dB)", legend=:bottomleft)
ps = []
for (i,(ri, ro, label, color, R)) in enumerate(zip(r_inner, r_outer, labels, colors, Rs))
    sheet = polyring(sides=4, s1=[D, 0], s2=[0, D], ntri=700, orient=45,
                     a=[ri], b=[ro], Rsheet=R, units=mm)
    push!(ps, plot(sheet, unitcell=true, title=label, lc=color))
    strata = [Layer()
              sheet
              Layer(width=5mm)
              pecsheet() # Perfectly conducting ground plane
              Layer()]
    results = analyze(strata, 1:0.2:25, (ϕ=0, θ=0), showprogress=false,
                      resultfile=devnull, logfile=devnull)
    data = extract_result(results, @outputs FGHz s11dB(h,h))
    plot!(p, data[:,1], data[:,2], label="PSSFSS "*label, lc=color)
    dat = readdlm("../src/assets/costa_2014_" * lowercase(label) * "_reflection.csv", ',')
    plot!(p, dat[:,1], dat[:,2], label="Costa "*label, ls=:dash, lc=color)
end
plot(ps..., layout=(1,3))